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In probability density function (PDF) methods a transport equation is solved numer- 
ically to compute the time and space dependent probability distribution of several flow 
variables in a turbulent flow. The joint PDF of the velocity components contains informa- 
tion on all one-point one-time statistics of the turbulent velocity field, including the mean, 
the Reynolds stresses and higher-order statistics. We developed a series of numerical al- 
gorithms to model the joint PDF of turbulent velocity, frequency and scalar compositions 
for high-Reynolds-number incompressible flows in complex geometries using unstructured 
grids. Advection, viscous diffusion and chemical reaction appear in closed form in the PDF 
formulation, thus require no closure hypotheses. The generalized Langevin model (GLM) 
is combined with an elliptic relaxation technique to represent the non-local effect of walls 
on the pressure redistribution and anisotropic dissipation of turbulent kinetic energy. The 
governing system of equations is solved fully in the Lagrangian framework employing a large 
number of particles representing a finite sample of all fluid particles. Eulerian statistics are 
extracted at gridpoints of the unstructured mesh. Compared to other particle-in-cell ap- 
proaches for the PDF equations, this methodology is non-hybrid, thus the computed fields 
remain fully consistent without requiring any specific treatment. Two testcases demon- 
strate the applicability of the algorithm: a fully developed turbulent channel flow and the 
classical cavity flow both with scalars released from concentrated sources. 

I. Introduction 

Probability density function (PDF) methods 1 ' 2 belong to the broader family of statistical approaches 
of turbulence modeling. As opposed to moment closure techniques, in PDF methods the full PDF of the 
turbulent flow variables is sought, which provides all one-point one-time statistical moments of the underlying 
fields. Raising the description to higher levels has several advantages. Convection and mean pressure appear 
in closed form and are treated mathematically exactly. The closure problem can be severe in combustion 
engineering, where developing accurate closure techniques for the highly nonlinear chemical source terms 
has proved elusive for realistic configurations. In large eddy simulation (LES), where most of the energy 
containing motions arc sought to be exactly resolved, at high Reynolds number and Damkohler number the 
chemical reactions take place at subgrid scales. Therefore these processes require closure assumptions in LES 
as well and the results have a first order dependence on the accuracy of these models. 3 In the PDF equations 
the source terms due to chemical reactions appear in closed form, thus no closure assumptions are necessary. 
Advection and viscous diffusion, processes that are fundamental in near-wall turbulent flows, are also in 
closed form. Closure hypotheses are necessary for the effect of fluctuating pressure, dissipation of turbulent 
kinetic energy and small-scale mixing of the scalar. In principle, a more complete statistical description is 
possible by solving for the full PDF instead of its specific moments as is done in moment closures. 

The price to pay for the increased level of description is the high dimensionality of the governing trans- 
port equation. As a consequence, instead of employing traditional numerical techniques, such as the finite 
difference or finite clement methods, the Eulerian field equations are written in a Lagrangian form and 
Monte Carlo methods are used to integrate a set of stochastic differential equations. Numerically, the flow 
is represented by a large number of Lagrangian particles that represent a finite sample of all fluid particles 
of the turbulent flow. This methodology not only has the advantage that Monte Carlo techniques are more 
economical for problems with high dimensionality, but the Lagrangian equations also appear in a significantly 
simpler form than their Eulerian counterparts. 
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A natural way of combining existing Eulcrian flow solvers with PDF methods is to develop hybrid 
formulations. Several authors reported on hybrid finite volume (FV)/Monte Carlo algorithms employing 
both structured 4 ' 5 and unstructured grids. 6,7 Different types of hybrid algorithms are possible depending 
on which quantities are computed in the Eulerian and Lagrangian framework and how the information 
exchange is carried out between the two representations. 4 However, a common characteristic of these hybrid 
formulations is that certain consistency conditions have to be met (and enforced on the numerical level) to 
ensure that all computed fields remain consistent throughout the simulation. 

We propose a non-hybrid formulation that is self-consistent both at the level of turbulence closure and the 
numerical method. Since no fields are computed redundantly, no consistency conditions have to be enforced. 
The joint PDF of velocity, characteristic turbulent frequency and a set of passive or reactive scalars are 
computed fully in the Lagrangian framework. An unstructured grid is used (i) to extract Eulerian statistics 
from particles, (ii) to solve for inherently Eulcrian quantities, such as the mean pressure and (Hi) to track 
particles throughout the flow domain. To model the high anisotropy and inhomogeneity of the Reynolds 
stress tensor in the vicinity of walls Drccben & Pope 8 combined Durbin's elliptic relaxation technique 9 
with the generalized Langevin model (GLM) of Haworth & Pope. 10 We have implemented the model in 
a general two-dimensional setting for complex geometries. The flow is resolved down to the viscous wall 
region, by imposing only the no-slip condition on particles without any damping or wall-functions. Two 
simple testcases are presented: a fully developed turbulent flow in a long-aspect-ratio channel geometry and 
a turbulent cavity flow, both with scalar releases. 



II. Governing equations 



T 



HE Eulcrian governing equation for a viscous incompressible flow is 

dt J dxj p dxi 



1^ + U 3 7^ + -— = ^ 2 U l , (1) 



where Ui, P, p and v are the Eulcrian velocity, pressure, constant density and kinematic viscosity, respectively. 
The Navier-Stokes equation (1) is written in the Lagrangian framework as a system of governing equations 
for Lagrangian particle locations Xi and velocities Ui 11 

dX, = Uidt+ (2v) 1/2 dW t (2) 

dZ*(f) = -^^dt + 2„ 7 ^dt+(2^ 2 ^dW J , (3) 

p OXi OXjOXj OXj 

where the isotropic Wiener process 12 dWi, which is a known stochastic process with zero mean and vari- 
ance dt, is identical in both equations and the Eulerian fields on the right hand side are evaluated at the 
particle locations Xi. The momentum of the particles governed by Equations (2) and (3) accounts for both 
advection and diffusion in physical space with a Gaussian probability distribution. After applying Reynolds 
decomposition to the Eulerian velocity Ui = (Ui) + Ui and pressure P = (P) + p we adopt the generalized 
Langevin model (GLM) 10 to model the appearing unclosed terms and obtain the stochastic model equation 
for the particle velocity increment 

+ dj {Uj - (Uj)) dt + (C e) 1/2 dWl, 

where Gij is a second-order tensor function, Co is a positive constant, s denotes the rate of dissipation of 
turbulent kinetic energy and dW[ is another Wiener process. In general, it is assumed that the tensor Gij 
is a function of the Reynolds stress tensor (muj), the dissipation rate e and the mean velocity gradients 
d (Ui) /dxj. The last two terms of Equation (4) jointly model the processes of pressure redistribution and 
dissipation of turbulent kinetic energy. The specification of Gy determines a particular local closure. Since 
minimal requirements on Gy and Co automatically guarantee realizablity, Gy can be specified so that the 
stochastic equation (4) yields equivalent statistics with any popular Reynolds stress closure at the level 
of second moments. 13 For near- wall flows Gy may also be specified by an elliptic equation based on the 
analogy that the pressure in incompressible flows is governed by a Poisson equation. This results in a non- 
locally determined Reynolds stress tensor, where the low-Reynolds-number effects of the wall are felt solely 
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through the boundary conditions of this elliptic equation. Close to the wall, the elliptic operator affects 
the solution, while far from the wall the solution blends into a local Reynolds stress model. Details on the 
elliptic relaxation technique are described by Durbin in the Reynolds stress framework and by Dreeben & 
Pope 8 in the PDF framework. 

Equation (4) can be closed by providing length or timescalc information for the turbulence. In traditional 
moment closures a model equation is solved for the turbulent kinetic energy dissipation rate 14 e or for 
the mean characteristic turbulent frequency 15 (uj) with the definition of the dissipation rate as e = k{uj), 
where k = | (uiUi) denotes the turbulent kinetic energy. In pure Lagrangian PDF methods, however, an 
alternative Lagrangian approach has been preferred for the characteristic particle frequency uj. A model for 
inhomogeneous flows has been developed by van Slooten & Pope, 16 whose simplest form is 

duo = -C 3 (uj) {uj - (uj)) dt - S u (uj) ujdt + (2C3C4 (uj) 2 cj) 1 2 dW, (5) 

where S u is a source/sink term for the mean turbulent frequency 

V 

Su = C^2 — Cui — , (6) 

£ 

where V = — (uiUj) d (Ui) /dxj is the production of turbulent kinetic energy, dW is a scalar-valued Wiener- 
process, while C3, C4, C^i and are model constants. 

In many practical engineering problems, such as combustion and atmospheric dispersion, the transport 
and dispersion of passive and reactive scalars in turbulence is of fundamental importance. A remarkable 
feature of PDF models is that since the source/sink terms in the advection-diffusion equations governing 
these scalar quantities appear in closed form, they can be represented mathematically exactly, without 
closure assumptions even in turbulent flows. The Eulerian governing equations for a set of reactive scalars 
(f> a , a = 1 . . . n 

^f + u- V0 Q = rv 2 (/) Q + s Q (0) (7) 

are written in the Lagrangian framework for the instantaneous particle compositions i\j a 

d^ a = rV 2 </> Q di + S a (ip)dt, (8) 

where the source terms S a (ip) are closed and closure is needed for the molecular diffusion term. For simplicity, 
the molecular diffusivity T is taken to be constant, uniform and the same for each composition. To model the 
molecular diffusion of the scalars we adopt the interaction by exchange with the conditional mean (IECM) 
model 

dVa = ~ (i>a - (<f> a \V))dt + S a &)dt, (9) 

where < m is a mixromixing timescale, while (cj) a \V) = (cp a \U(x,t) = V) denotes the expected value of the 
mean concentrations conditional on the velocity. For more on the theoretical details, see the reviews compiled 
by Pope 1 and Dopazo. 2 

In summary, the flow is represented by a large number of Lagrangian particles; their position Xi, velocity 
Ui, characteristic frequency u> and scalar concentrations ip a are governed by Equations (2), (4), (5) and 
(9), respectively. These equations are discretized and advanced in time with a forward Euler method. The 
Eulerian statistics appearing in the equations need to be evaluated at the particle positions. An Eulerian 
grid discretizes the flow domain and provides the spatial locations where these statistics are estimated. The 
mean pressure appearing in Equation (4) is computed by a pressure projection algorithm, also utilizing the 
Eulerian grid. 

III. Two testcases 

Two simple testcases demonstrate the applicability of the method: a fully developed turbulent channel 
flow and a more complex turbulent cavity flow. In the following two subsections selected statistics of 
the modeled joint PDF of velocity, frequency and a passive scalar are presented. 
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III.A. Channel flow 



The model has been run for a fully developed turbulent channel flow at the Reynolds number Re T = 1080 
based on the friction velocity u T and the channel half width h. After an initial development region this flow 
becomes statistically stationary, the velocity statistics become one-dimensional and remain inhomogeneous 
only in the wall-normal direction. Into this flow a passive scalar is released from a concentrated line source 
located at the channel centerline. Since the scalar field is inhomogeneous, the Eulerian grid is used to 
compute scalar statistics. Cross-stream profiles of (U), (uiUj) and e are plotted in Figure 1. Also shown are 
the DNS data of Abe et. al 17 at Re T = 1020. 

In turbulent channel flow the center region of the channel can be considered approximately homoge- 
neous. 18 Thus for a scalar released at the centreline, Taylor's theory of absolute dispersion 19 is expected to 
describe the mean field of the passive scalar well up to a certain downstream distance from the source. This 
is shown in Figure 2, where cross-stream mean concentration profiles at different downstream locations are 
depicted. Also shown in Figure 2 is a PDF of scalar concentration fluctuations (/>' = tp — (</>) at a location 
downstream of the source. The model for the joint PDF of U, u> and <\> accurately represents the full PDF 
and its statistics. 




y + 



Figure 1. Cross-stream profiles of (a) the mean streamwise velocity, (b) the diagonal components of the Reynolds 
stress tensor, (c) the shear Reynolds stress and (d) the rate of dissipation of turbulent kinetic energy. Lines — PDF 
calculation, symbols — DNS data of Abe et. al. 17 All quantities are normalized by the friction velocity and the channel 
half-width. The DNS data is scaled from Re T = 1020 to 1080. 



III.B. Cavity flow 

The turbulent cavity flow is used to demonstrate the applicability of the method in complex geometries. The 
model has been run at the Reynolds number Re = 4000 based on the free stream velocity and cavity depth. 
This flow also becomes statistically stationary after an initial period. After the flow is fully developed, a 
passive scalar is released from a concentrated source at the bottom of the cavity. The geometry of the domain 
and the spatial distributions of mean velocity, turbulent kinetic energy and the mean and variance of the 
scalar concentration are depicted in Figure 3. 
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Figure 2. (a) Cross-stream profiles of mean concentration at different downstream locations. Lines — PDF calculation, 
hollow symbols — analytical Gaussians according to Taylor's theory of dispersion, 19 filled symbols — experimental data 
of Lavertu & Midlarsky. 20 (b) PDF of concentration fluctuations at a downstream location at the centreline. Lines — 
computation, symbols — experimental data. 
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IV. Conclusion 



We developed and implemented a series of numerical methods to compute the joint PDF of turbulent 
velocity, frequency and scalar concentrations for high-Rcynolds-number incompressible turbulent flows with 
complex geometries. Adequate wall-treatment on the higher-order statistics is achieved with an elliptic 
relaxation technique without damping or wall functions. The current examples demonstrate the applicability 
of the algorithm for two dimensional flows. 
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